Hands-on Exercise 2

Author

PHAM Hung Son

Published

September 3, 2025

Modified

September 3, 2025

Part 1: First-order Spatial Point Patterns Analysis Methods

Importing and Wrangling Data

pacman::p_load(sf, terra, spatstat, 
               tmap, rvest, tidyverse)
mpsz_sf <- st_read("C:/sonphamsmu/isss626-aug2526/Hands-on_Ex/Hands-on_ex02/data/MasterPlan2019SubzoneBoundaryNoSeaKML.kml") %>% 
  st_zm(drop = TRUE, what = "ZM") %>% st_transform(crs = 3414)
Reading layer `URA_MP19_SUBZONE_NO_SEA_PL' from data source 
  `C:\sonphamsmu\isss626-aug2526\Hands-on_Ex\Hands-on_ex02\data\MasterPlan2019SubzoneBoundaryNoSeaKML.kml' 
  using driver `KML'
Simple feature collection with 332 features and 2 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: 103.6057 ymin: 1.158699 xmax: 104.0885 ymax: 1.470775
Geodetic CRS:  WGS 84
extract_kml_field <- function(html_text, field_name) {
  if (is.na(html_text) || html_text == "") return(NA_character_)
  
  page <- read_html(html_text)
  rows <- page %>% html_elements("tr")
  
  value <- rows %>%
    keep(~ html_text2(html_element(.x, "th")) == field_name) %>%
    html_element("td") %>%
    html_text2()
  
  if (length(value) == 0) NA_character_ else value
}
mpsz_sf <- mpsz_sf %>%
  mutate(
    REGION_N = map_chr(Description, extract_kml_field, "REGION_N"),
    PLN_AREA_N = map_chr(Description, extract_kml_field, "PLN_AREA_N"),
    SUBZONE_N = map_chr(Description, extract_kml_field, "SUBZONE_N"),
    SUBZONE_C = map_chr(Description, extract_kml_field, "SUBZONE_C")
  ) %>%
  select(-Name, -Description) %>%
  relocate(geometry, .after = last_col())
mpsz_cl <- mpsz_sf %>%
  filter(SUBZONE_N != "SOUTHERN GROUP",
         PLN_AREA_N != "WESTERN ISLANDS",
         PLN_AREA_N != "NORTH-EASTERN ISLANDS")
write_rds(mpsz_cl, 
          "C:/sonphamsmu/isss626-aug2526/Hands-on_Ex/Hands-on_ex02/data/mpsz_cl.rds")
childcare_sf <- st_read("C:/sonphamsmu/isss626-aug2526/Hands-on_Ex/Hands-on_ex02/data/ChildCareServices.kml") %>% 
  st_zm(drop = TRUE, what = "ZM") %>%
  st_transform(crs = 3414)
Reading layer `CHILDCARE' from data source 
  `C:\sonphamsmu\isss626-aug2526\Hands-on_Ex\Hands-on_ex02\data\ChildCareServices.kml' 
  using driver `KML'
Simple feature collection with 1925 features and 2 fields
Geometry type: POINT
Dimension:     XYZ
Bounding box:  xmin: 103.6878 ymin: 1.247759 xmax: 103.9897 ymax: 1.462134
z_range:       zmin: 0 zmax: 0
Geodetic CRS:  WGS 84
tm_shape(mpsz_cl) +
  tm_polygons() +
  tm_shape(childcare_sf) +
  tm_dots()

tmap_mode('view')
tm_shape(childcare_sf)+
  tm_dots()
tmap_mode('plot')
childcare_ppp <- as.ppp(childcare_sf)
summary(childcare_ppp)
Marked planar point pattern:  1925 points
Average intensity 2.417323e-06 points per square unit

Coordinates are given to 11 decimal places

Mark variables: Name, Description
Summary:
     Name           Description       
 Length:1925        Length:1925       
 Class :character   Class :character  
 Mode  :character   Mode  :character  

Window: rectangle = [11810.03, 45404.24] x [25596.33, 49300.88] units
                    (33590 x 23700 units)
Window area = 796335000 square units
sg_owin <- as.owin(mpsz_cl)
class(sg_owin)
[1] "owin"
plot(sg_owin)

childcareSG_ppp = childcare_ppp[sg_owin]
childcareSG_ppp
Marked planar point pattern: 1925 points
Mark variables: Name, Description 
window: polygonal boundary
enclosing rectangle: [2667.54, 55941.94] x [21448.47, 50256.33] units

Clark-Evan Test for Nearest Neighbour Analysis

The test hypotheses are:

Ho = The distribution of childcare services are randomly distributed.

H1= The distribution of childcare services are not randomly distributed.

The 95% confident interval will be used.

Without CSR method:

clarkevans.test(childcareSG_ppp,
                correction="none",
                clipregion="sg_owin",
                alternative=c("clustered"))

    Clark-Evans test
    No edge correction
    Z-test

data:  childcareSG_ppp
R = 0.53532, p-value < 2.2e-16
alternative hypothesis: clustered (R < 1)

Result: Since R = 0.53532 < 1 and p-value < 2.2e-16 (significant), we reject the null hypothesis.

Conclusion: The childcare centers are not randomly distributed but show statistically significant clustering.

With CSR method:

clarkevans.test(childcareSG_ppp,
                correction="none",
                clipregion="sg_owin",
                alternative=c("clustered"),
                method="MonteCarlo",
                nsim=99)

    Clark-Evans test
    No edge correction
    Monte Carlo test based on 99 simulations of CSR with fixed n

data:  childcareSG_ppp
R = 0.53532, p-value = 0.01
alternative hypothesis: clustered (R < 1)

Result: R < 1 and p = 0.01 (significant at the 5% level.)

Conclusion: Again Reject H₀. The distribution of childcare centers is statistically significantly clustered, not random.

Business Communication

Above analysis shows that childcare centers in Singapore are not randomly scattered across the city. Instead, they are statistically clustered, with some areas having a high density of centers while other neighborhoods have few or none.

This points to potential gaps in childcare accessibility that may require policy interventions or incentives to encourage new facilities in underserved regions.

For private operators, the clustering suggests high competition in certain areas but also highlights opportunities in neighborhoods where childcare services are lacking.

Kernel Density Estimation Method

kde_SG_diggle <- density(
  childcareSG_ppp,
  sigma=bw.diggle,
  edge=TRUE,
  kernel="gaussian") 
plot(kde_SG_diggle)

summary(kde_SG_diggle)
real-valued pixel image
128 x 128 pixel array (ny, nx)
enclosing rectangle: [2667.538, 55941.94] x [21448.47, 50256.33] units
dimensions of each pixel: 416 x 225.0614 units
Image is defined on a subset of the rectangular grid
Subset area = 669941961.12249 square units
Subset area fraction = 0.437
Pixel values (inside window):
    range = [-6.584123e-21, 3.063698e-05]
    integral = 1927.788
    mean = 2.877545e-06
bw <- bw.diggle(childcareSG_ppp)
bw
   sigma 
295.9712 
childcareSG_ppp_km <- rescale.ppp(
  childcareSG_ppp, 1000, "km")
kde_childcareSG_km <- density(childcareSG_ppp_km,
                              sigma=bw.diggle,
                              edge=TRUE,
                              kernel="gaussian")
plot(kde_childcareSG_km)

bw.CvL(childcareSG_ppp_km)
   sigma 
4.357209 
bw.scott(childcareSG_ppp_km)
 sigma.x  sigma.y 
2.159749 1.396455 
bw.ppl(childcareSG_ppp_km)
   sigma 
0.378997 
bw.diggle(childcareSG_ppp_km)
    sigma 
0.2959712 
kde_childcareSG.ppl <- density(childcareSG_ppp_km, 
                               sigma=bw.ppl, 
                               edge=TRUE,
                               kernel="gaussian")
plot(kde_childcareSG_km, main = "bw.diggle")

plot(kde_childcareSG.ppl, main = "bw.ppl")

Different kernel methods:

plot(density(childcareSG_ppp_km, 
             sigma=0.2959712, 
             edge=TRUE, 
             kernel="gaussian"), 
     main="Gaussian")

plot(density(childcareSG_ppp_km, 
             sigma=0.2959712, 
             edge=TRUE, 
             kernel="epanechnikov"), 
     main="Epanechnikov")

plot(density(childcareSG_ppp_km, 
             sigma=0.2959712, 
             edge=TRUE, 
             kernel="quartic"), 
     main="Quartic")

plot(density(childcareSG_ppp_km, 
             sigma=0.2959712, 
             edge=TRUE, 
             kernel="disc"), 
     main="Disc")

Fixed and Adaptive KDE

kde_childcareSG_fb <- density(childcareSG_ppp_km,
                              sigma=0.6, 
                              edge=TRUE,
                              kernel="gaussian")
plot(kde_childcareSG_fb)

kde_childcareSG_ab <- adaptive.density(
  childcareSG_ppp_km, 
  method="kernel")
plot(kde_childcareSG_ab)

plot(kde_childcareSG_fb, main = "Fixed bandwidth")

plot(kde_childcareSG_ab, main = "Adaptive bandwidth")

Plotting cartographic quality KDE map

kde_childcareSG_bw_terra <- rast(kde_childcareSG_km)
class(kde_childcareSG_bw_terra)
[1] "SpatRaster"
attr(,"package")
[1] "terra"
cat("\n")
kde_childcareSG_bw_terra
class       : SpatRaster 
size        : 128, 128, 1  (nrow, ncol, nlyr)
resolution  : 0.4162063, 0.2250614  (x, y)
extent      : 2.667538, 55.94194, 21.44847, 50.25633  (xmin, xmax, ymin, ymax)
coord. ref. :  
source(s)   : memory
name        :         lyr.1 
min value   : -5.824417e-15 
max value   :  3.063698e+01 
unit        :            km 
crs(kde_childcareSG_bw_terra) <- "EPSG:3414"
kde_childcareSG_bw_terra
class       : SpatRaster 
size        : 128, 128, 1  (nrow, ncol, nlyr)
resolution  : 0.4162063, 0.2250614  (x, y)
extent      : 2.667538, 55.94194, 21.44847, 50.25633  (xmin, xmax, ymin, ymax)
coord. ref. : SVY21 / Singapore TM (EPSG:3414) 
source(s)   : memory
name        :         lyr.1 
min value   : -5.824417e-15 
max value   :  3.063698e+01 
unit        :            km 
tm_shape(kde_childcareSG_bw_terra) + 
  tm_raster(col.scale = 
              tm_scale_continuous(
                values = "viridis"),
            col.legend = tm_legend(
            title = "Density values",
            title.size = 0.7,
            text.size = 0.7,
            bg.color = "white",
            bg.alpha = 0.7,
            position = tm_pos_in(
              "right", "bottom"),
            frame = TRUE)) +
  tm_graticules(labels.size = 0.7) +
  tm_compass() +
  tm_layout(scale = 1.0)

First Order SPPA at the Planning Subzone Level

pg <- mpsz_cl %>%
  filter(PLN_AREA_N == "PUNGGOL")
tm <- mpsz_cl %>%
  filter(PLN_AREA_N == "TAMPINES")
ck <- mpsz_cl %>%
  filter(PLN_AREA_N == "CHOA CHU KANG")
jw <- mpsz_cl %>%
  filter(PLN_AREA_N == "JURONG WEST")
plot(st_geometry(pg), main = "Punggol")

plot(st_geometry(tm), main = "Tampines")

plot(st_geometry(ck), main = "Choa Chu Kang")

plot(st_geometry(jw), main = "Jurong West")

pg_owin = as.owin(pg)
tm_owin = as.owin(tm)
ck_owin = as.owin(ck)
jw_owin = as.owin(jw)
childcare_pg_ppp = childcare_ppp[pg_owin]
childcare_tm_ppp = childcare_ppp[tm_owin]
childcare_ck_ppp = childcare_ppp[ck_owin]
childcare_jw_ppp = childcare_ppp[jw_owin]
childcare_pg_ppp.km = rescale.ppp(childcare_pg_ppp, 1000, "km")
childcare_tm_ppp.km = rescale.ppp(childcare_tm_ppp, 1000, "km")
childcare_ck_ppp.km = rescale.ppp(childcare_ck_ppp, 1000, "km")
childcare_jw_ppp.km = rescale.ppp(childcare_jw_ppp, 1000, "km")
plot(unmark(childcare_pg_ppp.km), 
  main="Punggol")

plot(unmark(childcare_tm_ppp.km), 
  main="Tampines")

plot(unmark(childcare_ck_ppp.km), 
  main="Choa Chu Kang")

plot(unmark(childcare_jw_ppp.km), 
  main="Jurong West")

Clark and Evans Test on Choa Chu Kang area:

clarkevans.test(childcare_ck_ppp,
                correction="none",
                clipregion=NULL,
                alternative=c("two.sided"),
                nsim=999)

    Clark-Evans test
    No edge correction
    Z-test

data:  childcare_ck_ppp
R = 0.84097, p-value = 0.008866
alternative hypothesis: two-sided

Tampines area:

clarkevans.test(childcare_tm_ppp,
                correction="none",
                clipregion=NULL,
                alternative=c("two.sided"),
                nsim=999)

    Clark-Evans test
    No edge correction
    Z-test

data:  childcare_tm_ppp
R = 0.66817, p-value = 6.58e-12
alternative hypothesis: two-sided
plot(density(childcare_pg_ppp.km, 
             sigma=bw.diggle, 
             edge=TRUE, 
             kernel="gaussian"),
     main="Punggol")

plot(density(childcare_tm_ppp.km, 
             sigma=bw.diggle, 
             edge=TRUE, 
             kernel="gaussian"),
     main="Tempines")

plot(density(childcare_ck_ppp.km, 
             sigma=bw.diggle, 
             edge=TRUE, 
             kernel="gaussian"),
     main="Choa Chu Kang")

plot(density(childcare_jw_ppp.km, 
             sigma=bw.diggle, 
             edge=TRUE, 
             kernel="gaussian"),
     main="Jurong West")

Part 2: 2nd Order Spatial Point Patterns Analysis Methods

set.seed(1234)
G_CK = Gest(childcare_ck_ppp, correction = "border") 
plot(G_CK, xlim=c(0,500))

Performing Complete Spatial Randomness Test

To confirm the observed spatial patterns above, a hypothesis test will be conducted. The hypothesis and test are as follows:

Ho = The distribution of childcare services at Choa Chu Kang are randomly distributed.

H1= The distribution of childcare services at Choa Chu Kang are not randomly distributed.

The null hypothesis will be rejected if p-value is smaller than alpha value of 0.001.

Monte Carlo test with G-fucntion

G_CK.csr <- envelope(childcare_ck_ppp, Gest, nsim = 999)
Generating 999 simulations of CSR  ...
1, 2, 3, ......10.........20.........30.........40.........50.........60..
.......70.........80.........90.........100.........110.........120.........130
.........140.........150.........160.........170.........180.........190........
.200.........210.........220.........230.........240.........250.........260......
...270.........280.........290.........300.........310.........320.........330....
.....340.........350.........360.........370.........380.........390.........400..
.......410.........420.........430.........440.........450.........460.........470
.........480.........490.........500.........510.........520.........530........
.540.........550.........560.........570.........580.........590.........600......
...610.........620.........630.........640.........650.........660.........670....
.....680.........690.........700.........710.........720.........730.........740..
.......750.........760.........770.........780.........790.........800.........810
.........820.........830.........840.........850.........860.........870........
.880.........890.........900.........910.........920.........930.........940......
...950.........960.........970.........980.........990........
999.

Done.
plot(G_CK.csr)

Tampines area:

G_tm = Gest(childcare_tm_ppp, correction = "best") 
plot(G_tm)

To confirm the observed spatial patterns above, a hypothesis test will be conducted. The hypothesis and test are as follows:

Ho = The distribution of childcare services at Tampines are randomly distributed.

H1= The distribution of childcare services at Tampines are not randomly distributed.

The null hypothesis will be rejected is p-value is smaller than alpha value of 0.001.

G_tm.csr <- envelope(childcare_tm_ppp, Gest, correction = "all", nsim = 999)
Generating 999 simulations of CSR  ...
1, 2, 3, ......10.........20.........30.........40.........50.........60..
.......70.........80.........90.........100.........110.........120.........130
.........140.........150.........160.........170.........180.........190........
.200.........210.........220.........230.........240.........250.........260......
...270.........280.........290.........300.........310.........320.........330....
.....340.........350.........360.........370.........380.........390.........400..
.......410.........420.........430.........440.........450.........460.........470
.........480.........490.........500.........510.........520.........530........
.540.........550.........560.........570.........580.........590.........600......
...610.........620.........630.........640.........650.........660.........670....
.....680.........690.........700.........710.........720.........730.........740..
.......750.........760.........770.........780.........790.........800.........810
.........820.........830.........840.........850.........860.........870........
.880.........890.........900.........910.........920.........930.........940......
...950.........960.........970.........980.........990........
999.

Done.
plot(G_tm.csr)

Analysing Spatial Point Process Using F-Function

CCK area:

F_CK = Fest(childcare_ck_ppp) 
plot(F_CK)

To confirm the observed spatial patterns above, a hypothesis test will be conducted. The hypothesis and test are as follows:

Ho = The distribution of childcare services at Choa Chu Kang are randomly distributed.

H1= The distribution of childcare services at Choa Chu Kang are not randomly distributed.

The null hypothesis will be rejected if p-value is smaller than alpha value of 0.001.

Monte Carlo test with F-fucntion

F_CK.csr <- envelope(childcare_ck_ppp, Fest, nsim = 999)
Generating 999 simulations of CSR  ...
1, 2, 3, ......10.........20.........30.........40.........50.........60..
.......70.........80.........90.........100.........110.........120.........130
.........140.........150.........160.........170.........180.........190........
.200.........210.........220.........230.........240.........250.........260......
...270.........280.........290.........300.........310.........320.........330....
.....340.........350.........360.........370.........380.........390.........400..
.......410.........420.........430.........440.........450.........460.........470
.........480.........490.........500.........510.........520.........530........
.540.........550.........560.........570.........580.........590.........600......
...610.........620.........630.........640.........650.........660.........670....
.....680.........690.........700.........710.........720.........730.........740..
.......750.........760.........770.........780.........790.........800.........810
.........820.........830.........840.........850.........860.........870........
.880.........890.........900.........910.........920.........930.........940......
...950.........960.........970.........980.........990........
999.

Done.
plot(F_CK.csr)

Tampines area:

F_tm = Fest(childcare_tm_ppp, correction = "best") 
plot(F_tm)

To confirm the observed spatial patterns above, a hypothesis test will be conducted. The hypothesis and test are as follows:

Ho = The distribution of childcare services at Tampines are randomly distributed.

H1= The distribution of childcare services at Tampines are not randomly distributed.

The null hypothesis will be rejected is p-value is smaller than alpha value of 0.001.

The code chunk below is used to perform the hypothesis testing.

F_tm.csr <- envelope(childcare_tm_ppp, Fest, correction = "all", nsim = 999)
Generating 999 simulations of CSR  ...
1, 2, 3, ......10.........20.........30.........40.........50.........60..
.......70.........80.........90.........100.........110.........120.........130
.........140.........150.........160.........170.........180.........190........
.200.........210.........220.........230.........240.........250.........260......
...270.........280.........290.........300.........310.........320.........330....
.....340.........350.........360.........370.........380.........390.........400..
.......410.........420.........430.........440.........450.........460.........470
.........480.........490.........500.........510.........520.........530........
.540.........550.........560.........570.........580.........590.........600......
...610.........620.........630.........640.........650.........660.........670....
.....680.........690.........700.........710.........720.........730.........740..
.......750.........760.........770.........780.........790.........800.........810
.........820.........830.........840.........850.........860.........870........
.880.........890.........900.........910.........920.........930.........940......
...950.........960.........970.........980.........990........
999.

Done.
plot(F_tm.csr)

Analysing Spatial Point Process Using K-Function

CCK area:

K_ck = Kest(childcare_ck_ppp, correction = "Ripley") 
plot(K_ck, . -r ~ r, ylab= "K(d)-r", xlab = "d(m)")

To confirm the observed spatial patterns above, a hypothesis test will be conducted. The hypothesis and test are as follows:

Ho = The distribution of childcare services at Choa Chu Kang are randomly distributed.

H1= The distribution of childcare services at Choa Chu Kang are not randomly distributed.

The null hypothesis will be rejected if p-value is smaller than alpha value of 0.001.

The code chunk below is used to perform the hypothesis testing.

K_ck.csr <- envelope(childcare_ck_ppp, Kest, nsim = 99, rank = 1, glocal=TRUE)
Generating 99 simulations of CSR  ...
1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20,
21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40,
41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60,
61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76, 77, 78, 79, 80,
81, 82, 83, 84, 85, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96, 97, 98, 
99.

Done.
plot(K_ck.csr, . - r ~ r, xlab="d", ylab="K(d)-r")

Tampines area:

K_tm = Kest(childcare_tm_ppp, correction = "Ripley") 
plot(K_tm, . -r ~ r,       ylab= "K(d)-r", xlab = "d(m)",       xlim=c(0,1000))

K_tm.csr <- envelope(childcare_tm_ppp, Kest, nsim = 99, rank = 1, glocal=TRUE)
Generating 99 simulations of CSR  ...
1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20,
21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40,
41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60,
61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76, 77, 78, 79, 80,
81, 82, 83, 84, 85, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96, 97, 98, 
99.

Done.
plot(K_tm.csr, . - r ~ r,       xlab="d", ylab="K(d)-r", xlim=c(0,500))

Analysing Spatial Point Process Using L-Function

CCK area:

L_ck = Lest(childcare_ck_ppp, correction = "Ripley") 
plot(L_ck, . -r ~ r,       ylab= "L(d)-r", xlab = "d(m)")

L_ck.csr <- envelope(childcare_ck_ppp, Lest, nsim = 99, rank = 1, glocal=TRUE)
Generating 99 simulations of CSR  ...
1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20,
21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40,
41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60,
61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76, 77, 78, 79, 80,
81, 82, 83, 84, 85, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96, 97, 98, 
99.

Done.
plot(L_ck.csr, . - r ~ r, xlab="d", ylab="L(d)-r")